%
% Trade From Space : Counterfactuals
% May 2025, andreas.moxnes@bi.no



%% Counterfactual I: Expand Panama Canal

clear
workdir = '';
path(workdir,path);

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%
t = 8;              % theta
b1 = 1.14;
delta = 1/t;
alpha = 1734;       % The value per dwt of container trade 

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Import data
%%%%%%%%%%%%%%%%%%%%%%%%%%

% Expenditure
import = csvread(strcat(workdir,'/expenditure.csv'));    
% Variables are A_acid traffic_dwt Emanu_port (manufacturing expenditure) E_port (total expenditure)
Y = import(:,4)*1000;             % Raw data is in 1000 USD, 2015 prices
YW = sum(Y);
N = length(Y);

% Traffic
% Each row is a departure port, each column is an arrival port 
import = csvread(strcat(workdir,'/trafficmatrix.csv'));         % (D x A)
portid = import(:,1);        % 1st column is list of port IDs
XiDwt = import(:,2:end);
Xi = XiDwt*alpha;

% Panama canal indicator
import = csvread(strcat(workdir,'/XPanamamatrix.csv'));         % (D x A)
XPanama = import(:,2:end);

% Port names and country
import = importdata(strcat(workdir,'/portname.csv'),',');
portnames = import.textdata(:,2);
cty = import.textdata(:,3);
portnames = join([portnames cty],', ');

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve GE
%%%%%%%%%%%%%%%%%%%%%%%%%%

% Change in t
k = b1^(-delta);
that = XPanama*k;
that(that==0)=1;


% Pi is Pi_hat^-theta, P is P_hat^-theta
[P Pi cc] = AIS_fixedpoint_hat(Y,Xi,that,t);

% Nominal wages
w = Pi.^(1/(t+1));

% Check if Y_hat_world = 1 (numeraire)
Yhat = Pi.^(1/(t+1));
disp('Check that change in world income=1'); disp(sum((Y/YW).*Yhat));

% Real wages
realinc = w./P.^(-1/t);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Traffic and trade (both in changes and levels)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Traffic flows, in changes
Xihat = that.^(-t).*repmat(P,1,N).*repmat(Pi',N,1);
Xi1 = Xi.*Xihat;    % Counterfactual level

% Trade : Baseline
D = diag(Y + (sum(Xi,1)'+ sum(Xi,2))/2);
Cx = inv(D-Xi);
X0 = Cx.*repmat(Y,1,N).*repmat(Y',N,1);

% Trade : Counterfactual
Yhat = Pi.^(1/(t+1));
Y1 = Y.*Yhat;
D = diag(Y1 + (sum(Xi1,1)'+ sum(Xi1,2))/2);
Cx = inv(D-Xi1);
X1 = Cx.*repmat(Y1,1,N).*repmat(Y1',N,1);
Xhat = X1./X0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Change in tau and exposure to Panama canal (pi)
% We do this for all possible Panama crossings kl
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Implied change in tau
tauhat = (Xhat.*repmat(Pi./Yhat,1,N).*repmat((P./Yhat)',N,1)).^(-1/t);

% % Export data
writematrix(Xhat,strcat(workdir,'/Xhat_v2.csv')) 
writematrix(tauhat,strcat(workdir,'/tauhat_v2.csv')) 
trade0 = X0-diag(diag(X0));
trade0 = trade0*YW;
writematrix(trade0,strcat(workdir,'/X0_v2.csv')); 
writematrix(XPanama,strcat(workdir,'/XPanama_v2.csv'));

[r c] = find(XPanama==1);
clear pimat;
for i=1:length(r)
  k=r(i); l=c(i);
  pimat(:,:,i) = (tauhat./(repmat(tauhat(:,k),1,N).*that(k,l).*repmat(tauhat(l,:),N,1))).^t;
end

% The main Panama crossing is Balboa-Manzanillo
BalboaIdx = find(strcmp(portnames, 'BALBOA*, PA'));
ManzanilloIdx = find(strcmp(portnames, 'MANZANILLO*, PA'));
rr1 = [r==BalboaIdx c==ManzanilloIdx];      % Balboa->Manzanillo
rr2 = [r==ManzanilloIdx c==BalboaIdx];      % Manzanillo->Balboa
idx1 = find(sum(rr1')'==2);
idx2 = find(sum(rr2')'==2);
pi1 = pimat(:,:,idx1);
pi2 = pimat(:,:,idx2);
pimax = max(pimat,[],3);


%% Counterfactual I : Results

close all
savefig=0;

disp('--------- Results --------- ');

% Figure : Top real income change, by port
figure(1);
top = 20;
C = [portnames,num2cell(realinc),num2cell(diag(Xhat))];
Csort = sortrows(C,2);
xx = (cell2mat(Csort(end-top+1:end,2))-1)*100;
tmp = Csort(end-top+1:end,1);
lab = categorical(tmp);
lab = reordercats(lab,tmp);
bar(lab,xx); ylabel('Real income, % change');
if savefig 
    saveas(gcf,strcat(workdir,'/graph/realincome.eps'),'eps');
    saveas(gcf,strcat(workdir,'/graph/realincome.pdf'),'pdf');
end

% Weighted average real income change
wch = sum(realinc.*(Y/YW));
disp('Welfare gains, %, weighted average'); disp((wch-1)*100);
disp('Welfare gains, in bill USD'); disp((wch-1)*YW/1000000000);

PAidx = find(strcmp(cty, 'PA'));
realincPA = sum(realinc(PAidx).*Y(PAidx))/sum(Y(PAidx));
PAbillgain = (realincPA-1)*sum(Y(PAidx))/1000000000;

worldbillgain = (wch-1)*YW/1000000000;
disp('Welfare gains for Panama, in bill USD'); disp(PAbillgain)
disp('Welfare gains for the world, in bill USD'); disp(worldbillgain);
disp('Ratio'); disp(worldbillgain/PAbillgain);  % Gains 3x bigger!


% Figure: Histogram of the change in tau and trade
figure(5);
subplot(1,2,1); hist((Xhat(:)-1)*100,50); xlabel('Trade, % change');
subplot(1,2,2); hist((tauhat(:)-1)*100,50); xlabel('\tau_{ij}, % change');
if savefig
  saveas(gcf,strcat(workdir,'/graph/tau_trade.eps'),'eps');
  saveas(gcf,strcat(workdir,'/graph/tau_trade.pdf'),'pdf');
end

% Change in tau and trade, %
idxdiag = eye(N);
disp('Unweighted average tau and X, % change');
disp(nanmean((tauhat(~idxdiag)-1)*100));
disp(nanmean((Xhat(~idxdiag)-1)*100));
disp('Median tau and X, % change');
disp(nanmedian((tauhat(~idxdiag)-1)*100));
disp(nanmedian((Xhat(~idxdiag)-1)*100));
disp('Max tau and X, % change');
disp(nanmax((tauhat(~idxdiag)-1)*100));
disp(nanmax((Xhat(~idxdiag)-1)*100));
disp('Min tau and X, % change');
disp(nanmin((tauhat(~idxdiag)-1)*100));
disp(nanmin((Xhat(~idxdiag)-1)*100));
disp('Increase in world trade, %');
disp((sum(X1(~idxdiag))./sum(X0(~idxdiag))-1)*100)


%% Counterfactual II: Drop links one at a time : Network model

clear
workdir = '';
path(workdir,path);

savefig=0;
top = 20;

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%
t = 40;
b1 = 1.14;
delta = 1/t;
alpha = 1734;       % The value per dwt of container trade 

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Import data
%%%%%%%%%%%%%%%%%%%%%%%%%%

% Expenditure
import = csvread(strcat(workdir,'/expenditure.csv'));    
% Variables are A_acid traffic_dwt Emanu_port (manufacturing expenditure) E_port (total expenditure)
Y = import(:,4)*1000;             % Raw data is in 1000 USD, 2015 prices
YW = sum(Y);
N = length(Y);

% Traffic
% Each row is a departure port, each column is an arrival port 
import = csvread(strcat(workdir,'/trafficmatrix.csv'));         % (D x A)
portid = import(:,1);        % 1st column is list of port IDs
XiDwt = import(:,2:end);
Xi = XiDwt*alpha;
%Xi = (Xi+Xi')/2;    % Make symmetric

% Panama canal indicator
import = csvread(strcat(workdir,'/XPanamamatrix.csv'));         % (D x A)
XPanama = import(:,2:end);

% Port names and country
import = importdata(strcat(workdir,'/portname.csv'),',');
portnames = import.textdata(:,2);
cty = import.textdata(:,3);
portnames = join([portnames cty],', ');

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve GE
%%%%%%%%%%%%%%%%%%%%%%%%%%

% Change in t

[r c] = find(Xi>0);

for i=1:length(r)
  disp('Link number'); disp(i);  
  that = ones(N,N);
  that(r(i),c(i)) = Inf;
  %that(c(i),r(i)) = Inf;       % Symmetric
  [P Pi cc] = AIS_fixedpoint_hat(Y,Xi,that,t);
  Pmat(:,i) = P;
  Pimat(:,i) = Pi;
  iter(i) = cc;
  % Pi is Pi_hat^-theta, P is P_hat^-theta
  
  w(:,i) = Pi.^(1/(t+1));               % Nominal wages
  realinc(:,i) = w(:,i)./P.^(-1/t);     % Real wages
  
  % Weighted average real income change
  wch(i) = sum(realinc(:,i).*(Y/YW));

end


writematrix(realinc,strcat(workdir,'/droplink_realinc_theta40.csv')) 
writematrix(Pmat,strcat(workdir,'/droplink_Pmat_theta40.csv')) 
writematrix(Pimat,strcat(workdir,'/droplink_Pimat_theta40.csv')) 
writematrix(wch,strcat(workdir,'/droplink_wch_theta40.csv')) 


%% Counterfactual II: View results

savefig=0;
top = 20;


realinc=readmatrix(strcat(workdir,'/droplink_realinc_theta40.csv')); 
Pmat=readmatrix(strcat(workdir,'/droplink_Pmat_theta40.csv')); 
Pimat=readmatrix(strcat(workdir,'/droplink_Pimat_theta40.csv')); 
wch=readmatrix(strcat(workdir,'/droplink_wch_theta40.csv')); 


[r c] = find(Xi>0);
if sum(isnan(realinc))>0 warning('Found NaN in realinc'); end

clear portpairs realinc_kl realinc_other
for i=1:length(r)
  portpairs(i,1) = join([portnames(r(i)) portnames(c(i))],' -> ');
end

% Real income change for sender and receiver (col 1 is sender, col 2 is
% receiver)
for i=1:length(r)
  realinc_kl(i,:) = [realinc(r(i),i) realinc(c(i),i)];
  realinc_kl_avg(i,1) = (realinc(r(i),i)*Y(r(i)) + realinc(c(i),i)*Y(c(i)))/(Y(r(i))+Y(c(i)));
  realinc_kl_bill(i,1) = (realinc_kl_avg(i,1)-1)*(Y(r(i))+Y(c(i)))/1000000000; % In bill
   
end

% Real income change for everybody else

for i=1:length(r)
  tmp = realinc(:,i);
  tmp([r(i) c(i)]) = [];
  Ytmp = Y;
  Ytmp([r(i) c(i)]) = []; 
  realinc_other(i,1) = sum(tmp.*(Ytmp/sum(Ytmp)));  
end

realinc_bill = sum((realinc-1).*repmat(Y,1,length(r)))'/1000000000;  % Change in world income, bill

% Total
C = [portpairs,num2cell((wch'-1)*100)];
Csort = sortrows(C,2,'descend');
xx = cell2mat(Csort(end-top+1:end,2));
tmp = Csort(end-top+1:end,1);
lab = categorical(tmp);
lab = reordercats(lab,tmp);


C = [portpairs,num2cell((realinc_other-1)*100)];
Csort = sortrows(C,2,'descend');
xx = cell2mat(Csort(end-top+1:end,2));
tmp = Csort(end-top+1:end,1);
lab = categorical(tmp);
lab = reordercats(lab,tmp);
figure(2); bar(lab,xx); ylabel('Real income, excluding kl, % change');

if savefig 
    saveas(gcf,strcat(workdir,'/graph/realincome_drop_other.eps'),'eps');
    saveas(gcf,strcat(workdir,'/graph/realincome_drop_other.pdf'),'pdf');
end


% Global change vs local change in real expenditure
relALL = realinc_bill./realinc_kl_bill;
disp('Share positive externalities'); disp(sum(relALL>1)/length(relALL));
disp('Average ratio'); disp(nanmean(relALL));
disp('Median ratio'); disp(nanmedian(relALL));

BalboaIdx = find(strcmp(portnames, 'BALBOA*, PA'));
ManzanilloIdx = find(strcmp(portnames, 'MANZANILLO*, PA'));
rr1 = [r==BalboaIdx c==ManzanilloIdx];      % Balboa->Manzanillo
rr2 = [r==ManzanilloIdx c==BalboaIdx];      % Manzanillo->Balboa
idx1 = find(sum(rr1')'==2);
idx2 = find(sum(rr2')'==2);

disp('Ratio Balboa->Manzanillo'); disp(relALL(idx1));  
disp('Ratio Manzanillo->Balboa'); disp(relALL(idx2));  

figure(3); 
relALL2 = relALL;
relALL2(relALL2>2) = 2;
relALL2(relALL2<.5) = .5;       % Censor extreme values
hist(relALL2,50);
set(gca, 'XLim', [.5, 2], 'XTick', .5:.5:2,'XTickLabel', .5:.5:2,'YLim', [0, 450]);

if savefig 
    saveas(gcf,strcat(workdir,'/graph/relloss_hist.eps'),'eps');
    saveas(gcf,strcat(workdir,'/graph/relloss_hist.pdf'),'pdf');
end


disp('% change global real income, excluding kl, mean'); disp(mean((realinc_other-1)*100));
disp('% change global real income, excluding kl, median'); disp(median((realinc_other-1)*100));
disp('Share negative growth'); disp(sum(realinc_other<1)./length(r));



%% Counterfactual III: Drop links one at a time : No-Network model

clear
workdir = '';
path(workdir,path);

savefig=0;
top = 20;

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%
t = 8;              % theta. If no heterogeneity (theta->Inf) then the algorithm converges to the shortest path
b1 = 1.14;
delta = 1/t;
alpha = 1734;       % The value per dwt of container trade 

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Import data
%%%%%%%%%%%%%%%%%%%%%%%%%%

% Expenditure
import = csvread(strcat(workdir,'/expenditure.csv'));    
% Variables are A_acid traffic_dwt Emanu_port (manufacturing expenditure) E_port (total expenditure)
Y = import(:,4)*1000;             % Raw data is in 1000 USD, 2015 prices
YW = sum(Y);
N = length(Y);

% Traffic
% Each row is a departure port, each column is an arrival port 
import = csvread(strcat(workdir,'/trafficmatrix.csv'));         % (D x A)
portid = import(:,1);        % 1st column is list of port IDs
XiDwt = import(:,2:end);
Xi = XiDwt*alpha;


% Convert traffic to trade
X0 = trafficToTrade(alpha,XiDwt,Y);
X0 = (X0+X0')/2;        % Make symmetric
trade = X0-diag(diag(X0));
disp('Aggregate import share'); disp(sum(trade(:))/sum(X0(:)));    


% Panama canal indicator
import = csvread(strcat(workdir,'/XPanamamatrix.csv'));         % (D x A)
XPanama = import(:,2:end);

% Port names and country
import = importdata(strcat(workdir,'/portname.csv'),',');
portnames = import.textdata(:,2);
cty = import.textdata(:,3);
portnames = join([portnames cty],', ');


%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve GE
%%%%%%%%%%%%%%%%%%%%%%%%%%

% Change in t

[r c] = find(Xi>0);

for i=1:length(r)
  disp('Link number'); disp(i);  
  that = ones(N,N);
  that(r(i),c(i)) = Inf;
  that(c(i),r(i)) = Inf;

  [P Pi cc] = AIS_fixedpoint_hat_nonetwork(X0,that,t);
  Pmat(:,i) = P;
  Pimat(:,i) = Pi;
  iter(i) = cc;
  % Pi is Pi_hat^-theta, P is P_hat^-theta
  
  w(:,i) = Pi.^(1/(t+1));               % Nominal wages
  realinc(:,i) = w(:,i)./P.^(-1/t);     % Real wages
  
  % Weighted average real income change
  wch(i) = sum(realinc(:,i).*(Y/YW));

  % World income (numeraire)
  Yhat = Pi.^(1/(t+1));
  YWhat(i) = sum((Y/YW).*Yhat);


end

% Export data
writematrix(realinc,strcat(workdir,'/droplink_nonet_symm_realinc.csv')) 
writematrix(Pmat,strcat(workdir,'/droplink_nonet_symm_Pmat.csv')) 
writematrix(Pimat,strcat(workdir,'/droplink_nonet_symm_Pimat.csv')) 
writematrix(wch,strcat(workdir,'/droplink_nonet_symm_wch.csv')) 
writematrix(YWhat,strcat(workdir,'/droplink_nonet_symm_YWhat.csv')) 

%% Counterfactual III: View results

savefig=0;
symm=1;

if symm==0 
  realinc=readmatrix(strcat(workdir,'/droplink_nonet_realinc.csv')); 
  Pmat=readmatrix(strcat(workdir,'/droplink_nonet_Pmat.csv')); 
  Pimat=readmatrix(strcat(workdir,'/droplink_nonet_Pimat.csv')); 
  wch=readmatrix(strcat(workdir,'/droplink_nonet_wch.csv')); 
  YWhat=readmatrix(strcat(workdir,'/droplink_nonet_YWhat.csv')); 
else
  realinc=readmatrix(strcat(workdir,'/droplink_nonet_symm_realinc.csv')); 
  Pmat=readmatrix(strcat(workdir,'/droplink_nonet_symm_Pmat.csv')); 
  Pimat=readmatrix(strcat(workdir,'/droplink_nonet_symm_Pimat.csv')); 
  wch=readmatrix(strcat(workdir,'/droplink_nonet_symm_wch.csv')); 
  YWhat=readmatrix(strcat(workdir,'/droplink_nonet_symm_YWhat.csv')); 
end


[r c] = find(Xi>0);

clear portpairs realinc_kl realinc_other
for i=1:length(r)
  portpairs(i,1) = join([portnames(r(i)) portnames(c(i))],' -> ');
end

% Real income change for sender and receiver (col 1 is sender, col 2 is
% receiver)
for i=1:length(r)
  realinc_kl(i,:) = [realinc(r(i),i) realinc(c(i),i)];
  realinc_kl_avg(i,1) = (realinc(r(i),i)*Y(r(i)) + realinc(c(i),i)*Y(c(i)))/(Y(r(i))+Y(c(i)));
  realinc_kl_bill(i,1) = (realinc_kl_avg(i,1)-1)*(Y(r(i))+Y(c(i)))/1000000000; % In bill   
end

% Real income change for everybody else

for i=1:length(r)
  tmp = realinc(:,i);
  tmp([r(i) c(i)]) = [];
  Ytmp = Y;
  Ytmp([r(i) c(i)]) = []; 
  realinc_other(i,1) = sum(tmp.*(Ytmp/sum(Ytmp)));  
end

realinc_bill = sum((realinc-1).*repmat(Y,1,length(r)))'/1000000000;  % Change in world income, bill

% Total
C = [portpairs,num2cell((wch'-1)*100)];
Csort = sortrows(C,2,'descend');
xx = cell2mat(Csort(end-top+1:end,2));
tmp = Csort(end-top+1:end,1);
lab = categorical(tmp);
lab = reordercats(lab,tmp);


% Others
C = [portpairs,num2cell((realinc_other-1)*100)];
Csort = sortrows(C,2,'descend');
xx = cell2mat(Csort(end-top+1:end,2));
tmp = Csort(end-top+1:end,1);
lab = categorical(tmp);
lab = reordercats(lab,tmp);

relALL = realinc_bill./realinc_kl_bill;
disp('Share positive externalities'); disp(sum(relALL>1)/length(relALL));
disp('Average ratio'); disp(mean(relALL));
disp('Median ratio'); disp(median(relALL));

BalboaIdx = find(strcmp(portnames, 'BALBOA*, PA'));
ManzanilloIdx = find(strcmp(portnames, 'MANZANILLO*, PA'));
rr1 = [r==BalboaIdx c==ManzanilloIdx];      % Balboa->Manzanillo
rr2 = [r==ManzanilloIdx c==BalboaIdx];      % Manzanillo->Balboa
idx1 = find(sum(rr1')'==2);
idx2 = find(sum(rr2')'==2);

disp('Ratio Balboa->Manzanillo'); disp(relALL(idx1)); 
disp('Ratio Manzanillo->Balboa'); disp(relALL(idx2));  

figure(3); 
relALL2 = relALL;
relALL2(relALL2>6) = [];
relALL2(relALL2<.2) = [];
hist(relALL,20);
set(gca, 'XLim', [.5, 2], 'XTick', .5:.5:2,'XTickLabel', .5:.5:2);
if savefig 
  if symm
    saveas(gcf,strcat(workdir,'/graph/relloss_nonet_symm_hist.eps'),'eps');
    saveas(gcf,strcat(workdir,'/graph/relloss_nonet_symm_hist.pdf'),'pdf');
  else
    saveas(gcf,strcat(workdir,'/graph/relloss_nonet_hist.eps'),'eps');
    saveas(gcf,strcat(workdir,'/graph/relloss_nonet_hist.pdf'),'pdf');
  end      
end





